Terrestrial evidence for ocean forcing of Heinrich events and subglacial hydrologic connectivity of the Laurentide Ice Sheet

During the last glacial period, the Laurentide Ice Sheet (LIS) underwent episodes of rapid iceberg discharge, recorded in ocean sediments as “Heinrich events” (HEs). Two competing models attempt to describe the stimulus for HEs via either internal ice sheet oscillations or external ocean-climate system forcing. We present a terrestrial record of HEs from the northeastern LIS that strongly supports ocean-climate forcing. Subglacial carbonate precipitates from Baffin Island record episodes of subglacial melting coincident with the three most recent HEs, resulting from acceleration of nearby marine-terminating ice streams. Synchronized ice stream acceleration over Baffin Island and Hudson Strait is inconsistent with internal ice sheet oscillations alone and indicates a shared ocean-climate stimulus to coordinate these different glaciological systems. Isotopic compositions of these precipitates record widespread subglacial groundwater connectivity beneath the LIS. Extensive basal melting and flushing of these aquifers during the last HE may have been a harbinger for terminal deglaciation.


INTRODUCTION
Over the course of the last 120 thousand years (ka), continental ice volumes and glacial climate conditions peaked during the last glacial maximum (LGM; ca. 26.5 to 19 ka ago) (1,2), followed by collapse of the Eurasian ice sheets and North American Laurentide Ice Sheet (LIS) and transition to an interglacial climate. Superimposed on these orbitally paced climate trends were millennial-scale oscillations in climate and ice dynamics linked to interactions among the iceocean-atmosphere systems. For instance, the so-called Dansgaard-Oeschger (D-O) events, originally observed in the oxygen isotope records of Greenland ice cores, reflect episodes of pronounced local atmospheric warming (3). D-O warm phases, or interstadials, are characterized by warm Greenland air temperatures, limited continental ice and freshwater flux into the North Atlantic, and strong North Atlantic deepwater (NADW) formation and Atlantic meridional overturning circulation (AMOC). D-O cold phases, or stadials, are instead characterized by cooler Greenland air temperatures, enhanced continental ice and freshwater flux to the North Atlantic, and dampened NADW formation and AMOC. Superimposed over some stadials, but not all, are Heinrich events (HEs): pulses of widespread iceberg discharge into the North Atlantic, predominantly from the LIS, that coincide with particularly fresh North Atlantic surface waters and especially weak NADW and AMOC (3,4). HEs are apparently linked to stadial conditions and occurred only during stadials and correlated climate proxies (5).
HEs were first identified as intervals of increased lithic grain abundances in deep-sea sediment cores from the North Atlantic that reflect deposition of ice-rafted detritus (IRD) dispersed widely by icebergs (6,7). They recur quasiperiodically every ca. 7 ka, and provenance studies have traced the bulk of HE-associated IRD to the Hudson Strait, implicating the Hudson Strait Ice Stream (HSIS; Fig. 1) as the predominant source of HE iceberg discharge, although other IRD sources do occur and are more prominent in some HE layers (4). Despite the unmistakable presence of HEs in ocean sedimentary archives and correlated climate proxies, the mechanisms initiating HSIS surging and iceberg swarms in the North Atlantic are still under debate. The quasiperiodicity and heterogeneous behavior of HEs (4) have motivated a longstanding effort to unravel their spatiotemporal complexity with a coherent mechanistic model. The onset of HEs depend on some combination of "internal" ice sheet dynamics and "external" ocean-climate forcings that enhanced ice streaming and iceberg production by the HSIS.
MacAyeal (8) proposed a model to describe HE cyclicity that is regulated entirely by internal ice sheet thermodynamics. Under this model, the HSIS begins in a state of reduced thickness and frozen basal conditions, resulting in slow ice discharge that causes the LIS to thicken. Eventually, ice thickness becomes sufficient to drive basal melting, which lubricates soft deformable sediments at the bed, causing surging of the HSIS and concomitant iceberg discharge to the North Atlantic, until the LIS thins to the initial condition. While this internal dynamics model emphasizes the importance of HSIS basal conditions and ice streaming in facilitating rapid release of icebergs into the North Atlantic, the external ocean-climate forcing model underscores the temporal coherence of HEs with fluctuations in the ocean-climate system. The observation that AMOC weakening not only overlaps HE-IRD layers but both precedes and follows them suggests that HE iceberg surging is a response to an oceanclimate forcing, rather than an internally paced ice sheet process (9). Marcott et al. (10) showed that reduced AMOC circulation, as during a stadial, led to convective mixing of low-latitude warm waters that propagated into the North Atlantic subsurface, where intermediatedepth warming maxima coincided with HEs. In their model (10), these warm waters came into contact with and melted ice shelves, which undermined their capacity to buttress ice streams and resulted in ice stream surging and rapid ice loss. Given the evidence for persistent sea-ice and open-water conditions in Baffin Bay (11) and the Labrador Sea (12) rather than extended ice shelf cover, Bassis et al. (13) emphasized the vulnerability of ice stream grounding lines to ocean warming. By modeling HSIS surging in response to grounding line melting modulated by glacial isostatic adjustment, they reproduced the differential tempos of HEs relative to D-O oscillations (13). These ocean-climate forced models are corroborated by recent work showing that North Atlantic IRD deposition began during peak subsurface warming (14).
HEs entailed the loss of large volumes of continental ice as meltwater and icebergs (4,15,16). Of particular importance is the final HE (H1, with progressively earlier events increasing in number), which coincided with the beginning of the last glacial termination. H1 entailed multiple phases of complex ice-ocean interactions spanning as much as 4 ka (17) and may have actively exacerbated terminal ice sheet collapse (18). During the time frame of H1 iceberg discharge and IRD deposition, the  234 U (calculated as  234 U = 1000×[( 234 U/ 238 U) − 1], where parentheses denote the activity ratio) of Atlantic surface waters rose by several per mille, indicating the influx of high- 234 U terrestrial waters, likely sourced from subglacial reservoirs (19). The precise subglacial source of this high- 234 U deglacial run-off is unclear, but its probable subglacial provenance temporally correlates ice sheet processes with H1 iceberg surging and terminal LIS collapse. A better understanding of the subglacial source of these distinct waters would clarify if these events are related and, if so, what glaciological mechanisms connected them. Like the North Atlantic IRD layers that characteristically record HEs, the Atlantic  234 U record, in this case from deep-sea corals, is a marine record linked to LIS perturbations during episodes of climate change. Connecting these marine archives to specific glaciological systems and processes is crucial for accurately reconstructing the glacial and deglacial histories of ice sheets.
To date, our current understanding of LIS dynamics over the course of HE and D-O cycles come from mathematical simulations, ocean sediment cores, and distal climate proxy records (5,13,20).  (73), paleo-ice streams draining the Foxe Dome (pale blue with dark borders) (33), subglacial carbonate precipitates measured in this study (diamonds), and groundwater-fed periglacial carbonate precipitates with  234 U o >500‰ (circles) (47,(58)(59)(60). Black boxes in (A) indicate extents of (B) to (D). (C) Reconstructed ice surface elevation and surrounding land and bathymetry (relative to modern sea level) at 20 ka ago (35). (D) Change in reconstructed ice thickness between 20 and 12.5 ka ago (35). evidences of LIS instability during HEs. On Baffin Island, however, the arid polar climate and persistence of cold-based ice caps have preserved subglacial calcite mineral precipitates near the margin of the modern-day Barnes Ice Cap (BIC) and the Rimrock Hills region (Fig. 1). Prior studies of these precipitates verified their subglacial origin and dated their formation approximately contemporaneous with the LGM (22,23). These aqueous precipitates require the presence of liquid water supersaturated with respect to calcite at a location that has previously been interpreted to be beneath cold-based ice (24). Acceleration of marine-terminating ice streams that drained the Foxe Dome of the LIS into Baffin Bay (Fig. 1) may have produced sufficient basal shear heating (25,26) to drive the requisite subglacial melting. Here, we report U-Th dates of multiple episodes of aqueous carbonate precipitation on Baffin Island and assess the provenance of the calcite-forming waters with stable and radiogenic isotope proxies. We propose that these subglacial precipitates record changes in the basal thermal regime due to HE-related ice streaming, supported by both the time scales of carbonate precipitation and the subglacial hydrologic systems in which they formed. By contextualizing these carbonate-forming waters relative to potential subglacial aquifers of the LIS, we substantiate a model of an expansive and broadly connected LIS subglacial groundwater system that was evacuated into marine basins by extensive basal melting during the early stages of deglaciation and H1.

Synchrony of HEs and subglacial calcite formation on Baffin Island
We calculated U-Th dates of six subglacial carbonate precipitates from the BIC margin (n = 5; figs. S3 to S6) and Rimrock Hills region (n = 1; fig. S7 and see Supplementary Text). These dates record intermittent episodes of aqueous precipitation of calcite spanning 31 to 18 ka ago (Table 1 and Fig. 2E). Five samples record single episodes of subglacial calcite formation between 18 and 26 ka ago. One sample (M09-B176R) records multiple episodes of subglacial calcite formation and erosion at the BIC margin: a basal layer precipitated at ca. 31 ka ago and an upper layer that records calcite formation between 25 and 23 ka ago, with a >4.6-ka disconformity separating the two layers (Table 1 and fig. S6). The most precise calcite formation ages (±<1 ka at 95% confidence intervals) identify three episodes of subglacial carbonate formation on central Baffin Island at ca. 31, 25 to 23, and ca. 18.5 ka ago. Less precise ages overlap at least one of these episodes within 95% confidence intervals ( Table 1). The time scales of carbonate formation we report are consistent with reliable U-Th dates reported by Refsnider et al. (22) for subglacial precipitates from the same localities (Table 1 and Supplementary Text). Figure 2 compares the ages of these subglacial calcite-forming events with several paleoclimate proxies. Precise U-Th ages align closely with HEs, as recorded by several proxies (Fig. 2): extended (HE) stadials in the NGRIP  18 O record (20), HE-associated positive excursions of the Hulu cave  18 O record (5), and enhanced surface freshening in the North Atlantic recorded by planktic  18 O (27,28). Baffin Island calcite precipitation also coincides with the carbonate IRD peaks of H2 and H3 in two North Atlantic marine sediment cores (10,14). Both of these cores also record North Atlantic subsurface seawater temperatures calculated from the Mg/Ca ratios of foraminifera: one (10) at intermediate depths (1 to 2 km) overlapping the intervals of H1 and H3 (H2 occurs during a gap in the record) and the other (14) at shallower subsurface depths (∼150 m) over the last 27 ka, overlapping H1 and H2. Baffin Island subglacial calcite formation coincides with subsurface warming during H1 and H3 and immediately follows the subsurface warming preceding H2 (Fig. 2, E, G, and H). Apparent lags between the onset of detrital carbonate deposition in the western North Atlantic recorded in marine sediment cores EW9302-2JPC (Fig. 2F) and GeoB18530-1 (Fig. 2I) and the Hulu Cave speleothem-constrained onset of H2 and H3 may reflect a lag between HSIS iceberg flux and global climate response or a systematic offset in radiocarbon and U-Th chronologies. Regardless, the time frames of subglacial calcite formation on Baffin Island overlap with H2 and H3, whether measured by U-Thconstrained global climate response or radiocarbon-constrained IRD deposition (Fig. 2).
The various records indicate concordance between Baffin Island carbonate precipitation events and HEs, implying a relationship between HE processes and basal ice sheet processes over central Baffin Island when it was covered by the eastern flank of the Foxe Dome of the LIS. Precipitation of calcite requires liquid water at or above the saturation point with respect to calcite. Yet, ice sheet models predict persistent cold-based conditions across Baffin Island during the last glacial period (24), corroborated by geologic evidence of relatively nonerosive basal conditions (29). Given this regional distribution of cold-based conditions, the occurrence of subglacial precipitates across hundreds of meters at sites separated by tens of kilometers (23) and evidence for hydrologic connectivity to groundwaters from the Canadian   (22), traced with a probability density plot (PDP) calculated using (76). Bells indicate uncertainty distributions of the means (vertical lines) and heights scale with precision (only relative heights matter, absolute values omitted). (F) Carbonate ice-rafted detritus (IRD) and (G) water temperatures calculated from benthic foraminifera Mg/Ca (±1 analytical uncertainty) at an intermediate-depth site in the Labrador Sea (core EW9302-2JPC) (10). (H) North Atlantic subsurface (∼150 m depth) water temperatures calculated from Mg/Ca ratios of sinistral N. pachyderma (envelope denotes 95% confidence interval) and (I) bulk sediment Ca/Sr, a proxy for carbonate IRD (core GeoB18530-1) (14). (J) Seawater  234 U (±2 uncertainty) in the upper (<1.5 km depth) low-latitude North Atlantic Ocean, reconstructed from deep-sea corals (19).
interior (see the following sections) implies episodes of widespread basal melting over northeastern Baffin Island during HEs. The preservation of delicate surface features on the carbonates (fig. S1) confirms cold-based ice at other times, including the late Holocene. Glacial striae on the bedrock from which samples were collected and linear features on several of the samples themselves (fig. S1) indicate basal ice sliding in the direction of fjords that were occupied by ice streams at the time (30).
Our record requires that subglacial melting occurred beneath the northern LIS synchronously with the HSIS acceleration and iceberg discharge canonically associated with HEs. This synchrony strongly favors a shared climatic triggering mechanism to coordinate the response and is not consistent with a stimulus of internally paced binge/purge oscillations of ice stream surging (8). This assertion is a consequence of the importance of ice streams in regulating ice flow from the dome. In a modern ice sheet, changes in ice stream behavior can reconfigure local mass balance (31). Even under the binge/purge model (8), HSIS acceleration was controlled by thickening of the ice stream, not the tributary domes. Investigations of surging glaciers, representing a well-studied example of glacial oscillators, show that such glaciers do not go through synchronous cycles, even if they are proximal, close in size, and subject to similar climatic conditions (32). Our precipitate samples formed in the upstream regions of outlet glaciers that drained the Foxe Dome eastward through Scott Inlet and Buchan Gulf into Baffin Bay. These glaciological settings differ distinctly from the HSIS, which drained several domes (15) through a ≥100-km-wide trough, well over twice the size of the 10-to 40-km-wide main and tributary troughs of the Scott and Buchan ice streams (30,33). The most parsimonious explanation for the fact that these glacial features experienced simultaneous ice acceleration during H1, H2, and H3 (Fig. 2) is that they reacted to a common repeat climate forcing rather than that they experienced synchronous oscillations unrelated to climate changes.
Alternatively, models relating HEs to subsurface ocean warming invoke HSIS surging in response to grounding line melting and retreat (13) or collapse of ice stream-buttressing ice shelves (10). Baffin Island ice streams that formerly occupied the Scott and Buchan Troughs, adjacent to the BIC, reached LGM grounding line extents >1000 m below modern-day sea level (30), the same intermediate water depths that record benthic warming correlated with HEs (10). We propose that basal melting and carbonate precipitation at the BIC margin occurred when Baffin Island ice streams surged in response to subsurface warming, accelerating inland ice velocities and elevating basal temperatures via shear heating (Fig. 3, Supplementary  Text, and fig. S8). We estimate that ice acceleration initiated at Baffin Bay grounding lines took 250 to 650 years to propagate to the site of BIC calcite formation (Supplementary Text), a lag within the reported age uncertainties of calcite formation (Table 1) and consistent with the ca. 500-year duration of IRD deposition observed during HEs (4). Although the present chronologies lack the temporal resolution to evaluate the timing of subglacial calcite formation relative to the onset of HEs, both the Baffin Bay sedimentary record and patterns of deglacial sea level response on Baffin Island support the interpretation that the ice stream responses recorded by subglacial calcite formation were a primary response to subsurface ocean warming, rather than a secondary response to ocean-climate consequences of HEs (Supplementary Text). The enhanced ice flow velocity must be connected to Baffin Island ice stream surging rather than HSIS surging given the relative stability of the Foxe Dome during widespread LIS thinning associated with HSIS activity: The Keewatin and Quebec-Labrador Domes thinned substantially during the early deglacial period (20 to 15 ka ago), whereas the Foxe Dome remained relatively unchanged (Fig. 1D). Moreover, the time scales required to propagate thinning across the Foxe Dome ice divide to the site of calcite formation are too long to explain the tight coincidence of subglacial calcite formation with HEs (Supplementary Text).
Currently, the Baffin Island subglacial precipitate chronology does not identify any subglacial carbonate precipitation on Baffin Island preceding H3. One potential explanation is that reduced thickness of the LIS during MIS 3 (34, 35) may have prevented Baffin Island ice streaming and basal melting during earlier HEs. Recuperated ice thickness at the BIC margin site by H3 (Fig. 2D) may have been necessary to sufficiently insulate the bed to accommodate melting during HE shear heating and/or to advance ice stream grounding lines to sufficient depth for exposure to subsurface warming events. Alternatively, limited preservation of pre-H3 carbonate precipitates may account for their absence in the present dataset. The disconformity in sample M09-B176R (Table 1 and  and might be identified by continued investigation of the BIC margin subglacial precipitate record. Conversely, despite ongoing subsurface ocean warming during H1, subglacial calcite precipitation stopped by ca. 18 ka ago, before the onset of H1 IRD deposition (Fig. 2, E, F, and I). This may be an artifact of a small number of samples that fail to capture the full duration of calcite precipitation. The distinct morphological character, concordance of formation, and summit locations of the two 18-to 19-ka Rimrock Hills samples imply that their formation mechanism may be decoupled from the BIC margin precipitates and sensitive only to basal melting following the LGM. In either case, the only post-LGM BIC margin sample (M09-B184R; Table 1) predates the bulk of H1 iceberg discharge (17). While further investigation may reveal additional calcites formed across the H1 interval, we currently interpret the abbreviated H1 calcite precipitation to imply a shift in LIS dynamics and associated basal conditions during the glacial termination that began at this time (36).

Subglacial fluid provenance
BIC margin and Rimrock Hills precipitates are restricted to the lee sides of bedrock undulations, bedrock fractures, and local depressions near summits, strongly supporting the conclusion that they form from the supersaturation of CaCO 3 , cryoconcentrated by the process of basal freezing in localized low-pressure zones (Fig. 3) (22,23). The extremely high U concentrations of the calcite (40 to >100 g/g; data file S1) imply formation from a U-rich fluid, further supporting the role of cryoconcentration in increasing the ionic strength of the calcite-forming waters. Comparable calcite U concentrations are precedented in speleothems from permafrost environments (37), where cryoconcentration of source waters likely also occurs.
Carbonate stable isotope compositions of  18 O and  13 C vary modestly (Fig. 4), consistent with the observations of Refsnider et al. (23). Our data reinforce their conclusions that  18 O reflects O isotopically fractionated during calcite precipitation from H 2 O in equilibrium with the overlying basal ice and  13 C reflects C isotopically fractionated from soil organic matter with minor contributions from bedrock calcite and atmospheric CO 2 . Over the course of >10 ka of intermittent calcite precipitation,  13 C does not change systematically with time, implying that this C source remained stable over this time frame (Fig. 4C). Similarly, 87 Sr/ 86 Sr does not vary systematically with time, except within sample M09-B176R, which records modest 87 Sr/ 86 Sr evolution during its growth (Fig. 4E) (Fig. 4), suggesting a shift toward isotopically lighter subglacial waters from melting of englacial ice (23) (Fig. 5). The heterogeneity in compositions distributed within a modeled three-component envelope requires >2 endmember waters or heterogeneity within "endmembers." One minor component may be Sr and U leached from silicate detritus during carbonate digestion procedures. Even with these complexities, the 87 Sr/ 86 Sr- 234 U data provide constraints on the approximate isotopic compositions and relative concentrations of Sr and U in two primary endmembers, arbitrarily named I and II (Fig. 5).
The highly radiogenic 87 Sr/ 86 Sr signature of endmember I (Fig. 5 (42). In addition, deep shield brines have higher total dissolved solids (TDS) than shallower groundwaters and converge on CaCl 2 compositions (43), both of which may explain the higher Sr content of endmember II than endmember I.
We are not aware of modern aqueous  234 U measurements on Baffin Island, but the elevated  234 U o compositions of both endmembers provide robust evidence for groundwater provenance (Fig. 5). Elevated aqueous  234 U is associated with the time-dependent ejection of 234 Th, which rapidly decays to 234 U, from sediment surfaces following -decay of 238 U. Aqueous 234 U enrichment in sediment/ rock pore space scales with increased residence time, lower porosity, and lower aqueous U concentration (44). Thus, while  234 U > 0 is typical among terrestrial waters, the largest enrichments are observed in groundwaters (45), with deep groundwater  234 U compositions commonly several thousand per mille (46). Waters with  234 U > 2000‰, as required by endmember II, are found in both bedrock aquifers (44,46) and groundwater-permafrost systems (37,47). We expect permafrost to have similar or elevated  234 U compared to liquid groundwaters of similar depth, as permafrost and slow-flowing subpermafrost groundwaters can experience protracted residence in contact with rock and sediment. Moreover, permafrost efficiently fractures rock (48), which increases surface area and enhances 234 U injection, a process that may have been augmented by subglacial hydraulic pressures from the LIS (49). Thus, we conclude that subglacial groundwater aquifers are the most probable source of the inferred high- 234 U waters recorded by the calcites from Baffin Island. Since deeper-residing groundwaters exhibit higher  234 U than their shallower counterparts by as much as 1000‰ or more (44,50), we identify the lower- 234 U endmember I as a relatively shallow groundwater and the higher- 234 U endmember II as a more deeply derived groundwater, consistent with the respective provenances inferred from Sr isotopes.
Stable isotope compositions do not vary systematically with 87 Sr/ 86 Sr (Fig. 4), indicating that mixing groundwater sources had little to no effect on O and C systematics. The most parsimonious explanation is that the local subglacial meltwater was initially well mixed with shallow Baffin Island groundwater (endmember I) and was the volumetrically dominant component of the calcite-forming waters. Sr and U reflect admixture with small volumes of high-TDS groundwaters (endmember II) that left isotopic compositions of H 2 O and DIC unperturbed. Even if larger volumes were incorporated, the effect on  13 C and  18 O would be minor, since deep saline and brine groundwaters have low DIC contents (51) and  18 O similar to or lighter than local meteoric water: Saline groundwaters in northern central Canada (52) approach the  18 O of BIC basal ice (22).
In summary, between 31 and 18 ka ago, Baffin Island shallow groundwaters were well mixed with LIS basal meltwater (in equilibrium with the basal ice) and had DIC predominantly derived from oxidized soil organic matter. These shallow groundwaters mixed with small volumes of saline or briny groundwaters distally sourced from deep bedrock aquifers toward the continental interior, altering cation compositions but having no measurable effect on the isotopic composition of H 2 O or DIC. Unlike the BIC margin site, the Rimrock Hills sites were isolated from distally sourced (endmember II) groundwaters, as evidenced by  234 U o <800‰ and inconsistency with the BIC margin mixing relationship (Fig. 5), perhaps due to their locations at local summits (22,39) and other regional factors that limited groundwater influx.

Subglacial groundwaters of the LIS
Precipitates at the modern-day BIC margin record the presence of waters from 31 to 18 ka ago that are compositionally consistent with high-salinity groundwaters found at depths of tens to hundreds of meters in crystalline basement rocks of the Canadian shield (42,44,51). The presence of these high-density shield brines at the base of the  LIS requires a powerful physical mechanism to transport them to the subglacial surface. We propose that this was achieved by glaciohydraulic reorganization of groundwater flow beneath the LIS. The overlying weight of ice sheets drove infiltration of basal meltwaters beneath the ice sheet interior, pressurizing subglacial aquifers and generating a steep hydraulic gradient between the ice sheet interior and margin (53,54). Fresh meltwaters were driven deep into the subsurface and mixed with high-TDS shield brines, as evidenced by geochemical observations of Canadian groundwaters at depths up to 1300 m (54,55) and numerical simulations of LIS groundwater recharge (56). As the water flowed along hydraulic gradients to the ice sheet margin, the flow paths shallowed and brought deeply sourced shield brines toward the surface. Simulations of sub-LIS groundwater flow predict upwelling of these brines extending ∼100 km interior from the northeastern LIS margin (56). Despite predicted permafrost conditions at the BIC margin over this time frame (24,56,57), aqueous mixing (Fig. 5) and calcite precipitation required intermittent subglacial permafrost melting and connectivity with deep groundwaters (Fig. 3), perhaps facilitated by fracture networks and freezing-point depression from the introduction of high-salinity brines. Moreover, the time-dependent evolution of 87 Sr/ 86 Sr in M09-B176R (Fig. 4) suggests that this hydrologic connectivity supported ongoing subglacial fluid transport over the course of basal melting events. In addition to the Baffin Island record, elevated  234 U compositions are found widely in marginal zones of the LGM LIS (Fig. 1A). Speleothems from the southern LIS margin and within 20 to 50 m of the modern surface record  234 U o ranging from 800 to 8000‰, with the majority >2000‰ over the last 250 ka (47). Late Pleistocene speleothems from caves near the southeastern LIS margin record  234 U o compositions in excess of 1000‰ (58). Cryogenically precipitated carbonates from LIS basal waters near the northwestern margin record  234 U o ∼800‰ at ca. 19 ka ago (59). Speleothems in the southern Canadian Rockies record  234 U o in excess of 1000‰ (60). The high- 234 U source waters may have been either sourced from deep groundwaters or evolved locally through protracted residence of porewater in permafrost environments (61). Permafrost provenance is less probable at low elevations since groundwater discharge during glacial terminations (57) would have interrupted the requisite porewater residence required for substantial in situ 234 U enrichment (61). Nonetheless, the groundwater-permafrost distinction is functionally unimportant: The mechanisms of 234 U-enrichment by -recoil injection are identical in bedrock aquifer, subglacial groundwater, and permafrost settings and all produce high- 234 U subglacial and periglacial groundwaters (44,61,62).
Collectively, subglacial precipitates on Baffin Island and marginal speleothems and cryogenic carbonates from around the LIS formed discontinuously in time, indicating that melting events at all of these marginal sites were typically transient or episodic. Although widespread subglacial hydrologic connectivity brought deep interior waters toward the LIS margins, the margins themselves were typically cold-based with robust marginal permafrost systems. These subsolidus marginal conditions not only restricted the growth of carbonate precipitates to critical intervals of climatic and glaciologic activity but also retained groundwaters within the ice sheet margins during glacial conditions.

Subglacial groundwaters and early deglacial ocean uranium chemistry
Records of high- 234 U shallow and surfacing groundwaters encircling the LIS margin ( Fig. 1) imply that such waters were commonplace, constituting a substantial reservoir of continental waters with elevated  234 U stored in permafrost and (potentially pressurized) subglacial aquifers (53). However, extensive subglacial melting could breach these confining cold-based margins and marginal permafrosts systems and efficiently route the high- 234 U waters to the margin. Although the areal extent of the LIS remained stable from the LGM until ca. 15 ka ago (1) and ice sheet models predict that LGM-like ice volumes persisted up to this time, areal basal melt increased from <40 to >60% between 25 and 15 ka ago (24), providing an efficient mechanism to melt permafrost, enhance hydrologic connectivity, and deliver interior high- 234 U waters to the ice sheet margin and ultimately the ocean.
Such a scenario of subglacial drainage would account for the abbreviated duration of H1 subglacial calcite formation. High- 234 U calcite-forming subglacial waters were present near the modern BIC margin during three consecutive HEs, but carbonate precipitation ended prematurely during H1, just as Atlantic  234 U compositions were rising toward their deglacial peak (Fig. 2, E and J). Before ca. 18 ka ago, we propose that calcite-forming waters were restrained by cold-based margins and continuous marginal permafrosts that persisted during nonterminal HEs. Then, during extended H1-deglacial conditions, enhanced LIS basal melting and hydrologic connectivity (24) breached these marginal features and routed the high- 234 U calcite-forming waters to the ice sheet margin, depleting the subglacial reservoir and shutting down carbonate growth by ca. 17 ka ago.
Therefore, we identify LIS groundwaters, broadly encompassing liquid porewaters in subglacial tills as well as permafrost in sub-and proglacial settings, as a compelling candidate for the subglacial reservoir responsible for the 6‰ enrichment of Atlantic Ocean  234 U during early deglaciation (19). Chen et al. (19) simulated the upper North Atlantic record by reducing a simplified model AMOC by 50% and enhancing the  234 U of modern riverine and groundwater inputs into the upper Atlantic from 339 to 800‰. In comparison, the  234 U of LIS-associated groundwaters considered here typically exceed 800‰. Subglacial precipitates do not reliably record the U concentration of their parent groundwaters; however, deep groundwaters with elevated  234 U are typically >10 g U/liter (46) and Arctic permafrosts with  234 U approaching 1000‰ range from ∼1 to 50 g/liter (61). Both exceed the U content of seawater (∼3 g/liter) (63), indicating that subglacial groundwaters and melted permafrost may serve as powerful levers on ocean  234 U at lesser volumes than riverine input (typically ≪3 g/liter) (45).
Marine and terrestrial records indicate that large-scale release of 234 U-rich subglacial fluids into the North Atlantic was apparently restricted to H1 and the early deglacial period. While H1 bore many of the hallmarks of preceding HEs, this ultimate iceberg discharge event was unique in both its evacuation of subglacial aquifers and its concurrence with termination of the last glacial period, implying a potential mechanistic relationship between these events (18,64). We speculate that before the LGM, marginal cold-based conditions and permafrost systems effectively restrained these interior subglacial aquifers. Following the LGM, extensive subglacial melting (24) enhanced hydrologic connectivity and degraded marginal permafrost systems, releasing subglacial waters that ultimately drained into the ocean. The concurrence of this process with H1 iceberg discharge implies that post-LGM subglacial melting promoted iceberg flux, HE-associated ice stream acceleration promoted extensive basal melting, or both processes reinforced each other. For now, the evidence is still limited and motivates further investigation of both subglacial and farther-afield records to examine the relationships among H1, the last glacial termination, and the role of subglacial hydrology in both.

Recontextualizing HEs
Subglacial calcites at the BIC margin formed during LIS basal melting events caused by nearby ice stream acceleration and implicate subsurface North Atlantic warming as a causal mechanism of HEs ( Fig. 3), independent of traditional HSIS evidences. Between HEs, the grounding lines of these ice stream systems readvanced to intermediate ocean water depths, requiring a rethickening of ice over Baffin Island, just as in the Hudson Strait. Although we show that ocean forcing initiated HEs, the recovery and regrowth phase, as in (8), remains a crucial component of the cycle. Notably, the coordinated ice stream response implies that Baffin ice streams and the HSIS recovered over similar time frames, potentially coupled via farreaching isostatic adjustment from eastern LIS domes (13) or rethickening of the shared Foxe Dome. Evidence for persistent open-water and sea-ice conditions over Baffin Bay (11) emphasizes the central role of ice streams, rather than extensive ice shelves, in controlling HE responses to subsurface ocean warming, though minor ice shelves and sea ice may have modulated these ice stream processes, e.g., via buttressing (10,65).
The climatic sensitivity of marine-terminating outlet glaciers on Baffin Island is supported by evidence of fluctuating sediment delivery from eastern Baffin Island into Baffin Bay on D-O time scales (66,67). While these sedimentary archives record outlet glacier responses to millennial-scale ocean-climate forcings, our data show that only HE conditions were sufficient to trigger an interior ice sheet response. Subglacial calcite forming events occurred on the interior highlands of Baffin Island, >600 m above modern sea level and >180 km from the LGM ice stream grounding lines. Since ice acceleration in response to grounding line retreat propagated at least this far into the ice sheet interior, this process reflects a regional, ice sheet-scale response to ocean forcing during HEs (Fig. 3). Yet, these ice sheet-scale episodes of acceleration and subglacial melting over Baffin Island were transient conditions: aqueous precipitates and nearby striated bedrock (23) require episodes of warm-based ice conditions, but inherited cosmogenic radionuclide signatures within the striated bedrock require a persistent state of cold-based, nonerosive conditions (68). To satisfy these observations, HE ice acceleration over Baffin Island must have reflected relatively brief (e.g., <1 ka) episodes of basal melting superimposed over the long-term cold-based conditions. The only sample that precisely records time across a single HE, the upper layer of sample M09-B176R (layers a1 to a5; Table 1), formed during H2 in as little as 800 years, corroborating the predicted brevity of melting episodes over the interior of Baffin Island. Substantial ice thicknesses over the region, ranging from ∼1.5 km at the site of calcite formation (Fig. 2D) to >3 km at the Foxe Dome (35), may have provided steep ice surface slopes that supported this rapid and extensive ice acceleration.
The rapid, coordinated regional response of ice streams terminating in both the Labrador Sea and Baffin Bay suggests that other maritime sectors of the eastern LIS, and perhaps other ice sheets, may have responded similarly to ocean forcing during HEs. Widespread subsurface ocean warming and subsequent surging of various marineterminating ice streams in the North Atlantic would account for the heterogeneity in IRD provenance observed in HE layers, which suggests iceberg flux from a variety of sources including Baffin Island, Greenland, Iceland, and Eurasian ice sheets (4,69). Such oceancoordinated ice sheet responses across the North Atlantic provide a framework to interpret the relationship between HEs and their northern Pacific analog, Siku events-massive iceberg discharge episodes of the Cordilleran Ice Sheet that precede HEs (70). On the basis of our findings of coordinated ice stream acceleration during HEs, we speculate that episodes of ice sheet-wide ice stream activation reflect intraoceanic responses to millenial-scale interoceanic climate oscillations.

MATERIALS AND METHODS
We studied six subglacially formed detritus-rich carbonate precipitate rocks, including five samples from the northern BIC margin and one sample from the Rimrock Hills region ( Fig. 1 and fig. S1), originally reported by Refsnider et al. (22,23). Small "crust" samples (M09-B184R and M09-B152) were subsampled with steel hand tools, while larger samples were slabbed and subsampled with a rock saw before isolating individual fractions with hand tools.

Carbon and oxygen isotope measurements
Isotopes of oxygen and carbon in carbonate phases were measured at the UC Santa Cruz Stable Isotopes Laboratory by acid digestion using an individual vial acid drop Thermo Fisher Scientific Kiel IV carbonate device interfaced to a Thermo Fisher Scientific MAT 253 dual-inlet isotope ratio mass spectrometer (iRMS). We loaded ∼100-g fractions into individual vials that were dried overnight in a 70°C vacuum oven. Samples reacted at 75°C with H 3 PO 4 (1.92 g/cm 3 specific gravity) to generate CO 2 and H 2 O. The latter is cryogenically separated, and noncondensible gases are pumped away before introduction of the CO 2 analyte into the iRMS. Uranium, thorium, and strontium isotope measurements Individual carbonate precipitate fractions were cleaned by sonication at room temperature in methanol for 30 min, triple-rinsed with ultrapure water (deionized to 18 megohms·cm), and transferred to an acidcleaned perfluoroalkoxy alkane (PFA) beaker. All aqueous reagents (except ultrapure water) were either triple-distilled or commercial tracemetal grade. We submerged cleaned samples in approximately 1 ml of ultrapure water and added 7 M HNO 3 dropwise to gently digest the carbonate minerals. When the reaction subsided, we brought the solution to 3 M HNO 3 , spiked with a gravimetrically calibrated mixed 229 Th-236 U tracer, and warmed the solution for several hours to progress reactions to completion. We twice evaporated the fractions to dryness, rehydrated in 7 M HNO 3 , and refluxed at 110°C to ensure complete digestion and sample-spike equilibration.
We purified Th and U separates from each dissolved fraction using ion chromatography on a 1-ml bed of AG1-X8 anion-exchange resin (200 to 400 mesh). We introduced dissolved samples onto pre-cleaned resin in 1 ml of 7 M HNO 3 and eluted matrix ions (including Sr) with 2 ml of 7 M HNO 3 followed by 250 l of 6 M HCl. We eluted Th in 2 ml of 6 M HCl followed by U in 3 ml of water into a single beaker. We twice evaporated the U-Th separate to dryness and rehydrated with 7 M HNO 3 to convert to nitrate salts. We then repeated the same ion chromatography procedure and collected the Th and U elutions separately. We evaporated the elutions just to dryness, refluxed at 110°C in 250 l of 30% H 2 O 2 for several hours to remove organic compounds, and lastly dried the solutions with trace H 3 PO 4 . Total procedural blanks of Th and U were <50 and <80 pg, respectively, and negligible compared to sample sizes (hundreds of nanograms of U and ≫10 ng of Th).
To purify Sr, we evaporated matrix elutions from the primary U-Th column to dryness, rehydrated in 0.5 ml of 7 M HNO 3 , and introduced this solution onto a 0.5-ml bed of precleaned Sr-Spec resin. We eluted matrix with 3 ml of 7 M HNO 3 and collected purified Sr in 4 ml of 0.05 M HNO 3 . We dried the latter with trace H 3 PO 4 .
We measured isotopes of U, Th, and Sr on the Isotopx X62 thermal ionization mass spectrometer at UC Santa Cruz. All samples were loaded onto degassed 99.99% purity Re ribbons. We loaded U with a Si gel-0.035 M H 3 PO 4 activator and measured UO 2 isotopologs using a dynamic Faraday-Daly method (62). We calculated U isotope compositions from UO 2 isotopologs by correcting for oxide isobaric interferences, although these corrections were negligible compared to analytical uncertainties. We corrected for mass-dependent fractionation with a linear model calibrated from long-term standard measurements and calibrated the photomultiplier deadtime from measurements of NBS SRM U-500 (as UO 2 ). The accuracy of U isotope measurements were validated over the course of this study with replicate measurements of National Institute of Standards and Technology Standard Reference Material (NIST SRM) 4321b ( fig. S2). We loaded Th with 1 l of 5% HNO 3 onto Re filaments coated with graphite and measured isotopes of Th as a metal on the Dalyphotomultiplier complex using a peak hopping method. Th isotope ratios were corrected for mass-dependent fractionation and photomultiplier deadtime using model values determined by measurements of SRM U-500 ionized as a metal.
U-Th isotope measurements were spike-subtracted with an algorithm that fully propagates analytical and tracer uncertainties, assuming uncorrelated uncertainties. We report and interpret measured isotopic data as activity ratios, denoted with parentheses. In the case of U isotopes, we report carbonate U compositions as activity ratios and report inferred calcite-forming water compositions in  notation. To precisely calculate U-Th dates and  234 U o , we used a Monte Carlo method algorithm (10 6 trials per fraction) that fully propagates the uncertainties of all input isotopic ratios, including corrections for detrital contributions where necessary. Dates are calculated relative to a 1950 CE "present" datum. The accuracy of U-Th dates and  234 U o was confirmed with concurrent measurement of a Marine Isotope Stage 5e coral (data file S1).
We loaded Sr with a TaCl 5 activator and measured Sr isotopes on Faraday cups with a static collection method. We correct for isobaric interference from 87 Rb on 87 Sr by concurrent measurement of 85 Rb and subtracting the intensity scaled by an assumed 87 Rb/ 85 Rb = 0.386. We calculate a mass-dependent fractionation  factor from the measured 86 Sr/ 88 Sr with an exponential law, assuming a canonical 86 Sr/ 88 Sr = 0.11940. We calculate fully corrected radiogenic 87 Sr/ 86 Sr compositions by applying this  factor to 87 Rb-corrected 87 Sr/ 86 Sr ratios. Before data collection, we confirmed reproduction of NIST SRM 987: 87 Sr/ 86 Sr = 0.710250 ± 0.000011 (1 SD) (72). Over the course of this study, we report a mean SRM 987 87 Sr/ 86 Sr = 0.710232 ± 0.000027 (2 SE, n = 8), consistent with the accepted value. Given the detritus-rich nature of the carbonate rock samples, we consider the U-Th-Sr isotope systematics in detail and quantify radiogenic and detrital components with isochron-based methods to calculate U-Th dates (Supplementary Text and figs. S3 to S7).